#include "routine_numeric.h"
/***************************
determinant of a real matrix
A: n x n
****************************/
double det(double *A, int n)
{
  double val=0; 
  double A0i=0;
  double *co_A=(double *)malloc(sizeof(double));
  co_A[0]=0; free(co_A);

  int i, ii, jj, ix, iy; 

  if(n==1){
    val=A[0];
  }
  else if(n==2){
    val=A[0+2*0]*A[1+2*1]-A[1+2*0]*A[0+2*1];
  }
  else{
    co_A=(double *)malloc((n-1)*(n-1)*sizeof(double)); 
    for(i=0;i<n;i++){
       A0i=A[i+n*0];

       for(ii=0;ii<n-1;ii++){
          for(jj=0;jj<n-1;jj++){

             if(ii<i){
               ix=ii;
             }
             else{
               ix=ii+1;
             }
               
             iy=jj+1;
             
             co_A[ii+(n-1)*jj]=A[ix+n*iy];
          }
       } 
       val= val+ pow(-1,i)*A0i*det(co_A,n-1); 
    } 
    free(co_A);
  }  
  return val;
}


/**
 * calculate the integral of a scalar grid function
 * Assuming that the grid space are uniform in any direction
 * the value of the grid_fun[iy][ix][iz]
 */
double complex Intl_grid_fun_1(double complex *grid_fun,
			       int Nx, int Ny, int Nz,
			       double hx, double hy, double hz)
{
  int ix, iy, iz;
  double complex Intl_val;  // the return value
  double complex u[8];
  Intl_val = 0.0 + 0.0*csqrt(-1.0);

  // 1. add the eight corner
  for(ix=0;ix<Nx-1;ix++){
    for(iy=0;iy<Ny-1;iy++){
      for(iz=0;iz<Nz-1;iz++){
	u[0] = grid_fun[iz+Nz*(ix+Nx*iy)];
	u[1] = grid_fun[iz+Nz*(ix+Nx*(iy+1))];
	u[2] = grid_fun[iz+Nz*(ix+1+Nx*iy)];
	u[3] = grid_fun[iz+Nz*(ix+1+Nx*(iy+1))];
	u[4] = grid_fun[iz+1+Nz*(ix+Nx*iy)];
	u[5] = grid_fun[iz+1+Nz*(ix+Nx*(iy+1))];
	u[6] = grid_fun[iz+1+Nz*(ix+1+Nx*iy)];
	u[7] = grid_fun[iz+1+Nz*(ix+1+Nx*(iy+1))]; 
	Intl_val += 0.125*(u[0] + u[1] + u[2] + u[3] + 
			   u[4] + u[5] + u[6] + u[7]);
      }
    }
  }
  Intl_val *= hx*hy*hz;
  return Intl_val;
}

// with coef, the integral is equation to inner product with coef
void Intl_grid_coef(double *coef, int Nx, int Ny, int Nz)
{
  int k, i, ix, iy, iz;
  // fill the first layer 2D
  k = 0;
  coef[k++] = 1.0;
  for(iz = 1; iz < Nz-1; iz++) coef[k++] = 2.0; 
  coef[k++] = 1.0;
  for(ix = 1; ix < Nx-1; ix++){
    coef[k++] = 2.0;
    for(iz = 1; iz < Nz-1; iz++)coef[k++] = 4.0;
    coef[k++] = 2.0;
  }
  coef[k++] = 1.0;
  for(iz = 1; iz < Nz-1; iz++) coef[k++] = 2.0; 
  coef[k++] = 1.0;

  // from the second layer to last second layer
  for(iy = 1; iy < Ny-1; iy++){
    coef[k++] = 2.0;
    for(iz = 1; iz < Nz-1; iz++) coef[k++] = 4.0; 
    coef[k++] = 2.0;
    for(ix = 1; ix < Nx-1; ix++){
      coef[k++] = 4.0;
      for(iz = 1; iz < Nz-1; iz++)coef[k++] = 8.0;
      coef[k++] = 4.0;
    }
    coef[k++] = 2.0;
    for(iz = 1; iz < Nz-1; iz++) coef[k++] = 4.0; 
    coef[k++] = 2.0;
  }

  // the last layer
  coef[k++] = 1.0;
  for(iz = 1; iz < Nz-1; iz++) coef[k++] = 2.0; 
  coef[k++] = 1.0;
  for(ix = 1; ix < Nx-1; ix++){
    coef[k++] = 2.0;
    for(iz = 1; iz < Nz-1; iz++)coef[k++] = 4.0;
    coef[k++] = 2.0;
  }
  coef[k++] = 1.0;
  for(iz = 1; iz < Nz-1; iz++) coef[k++] = 2.0; 
  coef[k++] = 1.0;

}


/**
 *   1. FFTW Code for -\laplacian u= 4*pi*rho
 *   
 */
void poissonFFT(double complex *v_H, double complex *rho, 
		int Nx, int Ny, int Nz, double hx, double hy, double hz)

{
  int ix=0, iy=0, iz=0;
  const int n0=Ny-1, n1=Nx-1, n2=Nz-1;
  int NNy=n0, NNx=n1, NNz=n2;
  double theta_y=2.0*PI/NNy, theta_x=2.0*PI/NNx, theta_z=2.0*PI/NNz;
  fftw_complex in[n0][n1][n2],out[n0][n1][n2];
  fftw_plan p;

  /****
       0. ready data
   ****/
  for(iy=0;iy<n0;iy++){
    for(ix=0;ix<n1;ix++){
      for(iz=0;iz<n2;iz++){
	in[iy][ix][iz]=4.0*PI*(rho[iz+Nz*(ix+Nx*iy)]);
	//in[iy][ix][iz][0]=4.0*PI*creal(rho[iz+Nz*(ix+Nx*iy)]);
	//in[iy][ix][iz][1]=4.0*PI*cimag(rho[iz+Nz*(ix+Nx*iy)]); 
	out[iy][ix][iz]=0;
      }
    }
  }

  /****
       1. forward FFT
   ****/
  p=fftw_plan_dft_3d(n0,n1,n2,&in[0][0][0],&out[0][0][0],FFTW_FORWARD,FFTW_ESTIMATE);
  fftw_execute(p); 
  fftw_destroy_plan(p); 

  /****
       2. Spectral Accuracy solving
   ****/  
  double denomy[n0];
  double denomx[n1];
  double denomz[n2];
  double denom=0;
  for(iy=0;iy<n0;iy++){
    if(iy<=NNy/2-1){
      denomy[iy]=iy;
    }
    else{
      denomy[iy]=iy-NNy;
    }
  } 
  for(ix=0;ix<n1;ix++){
    if(ix<=NNx/2-1){
      denomx[ix]=ix;
    }
    else{
      denomx[ix]=ix-NNx;
    }
  } 
  for(iz=0;iz<n2;iz++){
    if(iz<=NNz/2-1){
      denomz[iz]=iz;
    }
    else{
      denomz[iz]=iz-NNz;
    }
  } 

  
  for(iy=0;iy<n0;iy++){
    denomy[iy]=2.0*PI/(hy*NNy)*denomy[iy];
    denomy[iy]=denomy[iy]*denomy[iy];
  }
  for(ix=0;ix<n1;ix++){
    denomx[ix]=2.0*PI/(hx*NNx)*denomx[ix];
    denomx[ix]=denomx[ix]*denomx[ix];
  }
  for(iz=0;iz<n2;iz++){
    denomz[iz]=2.0*PI/(hz*NNz)*denomz[iz];
    denomz[iz]=denomz[iz]*denomz[iz];
  }
  
  for(iy=0;iy<n0;iy++){
    for(ix=0;ix<n1;ix++){
      for(iz=0;iz<n2;iz++){
	denom=denomy[iy]+denomx[ix]+denomz[iz];
	if(denom!=0){
	  out[iy][ix][iz]=out[iy][ix][iz]/denom;
	  //out[iy][ix][iz][0]=out[iy][ix][iz][0]/denom;
	  //out[iy][ix][iz][1]=out[iy][ix][iz][1]/denom;
	}
      }
    }
  }//*/ 
 
  /****
       3. IFFT u
   ****/
  p=fftw_plan_dft_3d(n0,n1,n2,&out[0][0][0],&in[0][0][0],FFTW_BACKWARD,FFTW_ESTIMATE); 
  fftw_execute(p);  
  fftw_destroy_plan(p); 

 
  /****
       4. return value
   ****/
  for(ix=0;ix<n1;ix++){
    for(iy=0;iy<n0;iy++){
      for(iz=0;iz<n2;iz++){ 
	v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][iz]/(NNx*NNy*NNz);  
      }
    }
  } 

  
  /****
       5. keep periodic?
  ****/
  for(ix=0;ix<Nx;ix++){
    for(iy=0;iy<Ny;iy++){
      for(iz=0;iz<Nz;iz++){
	if(ix==Nx-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=v_H[iz+Nz*(0+Nx*iy)];
	  //v_H[iz+Nz*(ix+Nx*iy)]=in[iy][0][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][0][iz][1]/(NNx*NNy*NNz);
	}
	else if(iy==Ny-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=v_H[iz+Nz*(ix+Nx*0)];
	  //v_H[iz+Nz*(ix+Nx*iy)]=in[0][ix][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[0][ix][iz][1]/(NNx*NNy*NNz);
	}
	else if(iz==Nz-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=v_H[0+Nz*(ix+Nx*iy)];
	  //v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][0][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][ix][0][1]/(NNx*NNy*NNz);
	}
	/*
	  if(ix<Nx-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][ix][iz][1]/(NNx*NNy*NNz);
	  }
	  else{
	  v_H[iz+Nz*(ix+Nx*iy)]=in[iy][0][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][0][iz][1]/(NNx*NNy*NNz);
	  }

	  if(iy<Ny-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][ix][iz][1]/(NNx*NNy*NNz);
	  }
	  else{
	  v_H[iz+Nz*(ix+Nx*iy)]=in[0][ix][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[0][ix][iz][1]/(NNx*NNy*NNz);
	  }

	  if(iz<Nz-1){
	  v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][iz][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][ix][iz][1]/(NNx*NNy*NNz);
	  }
	  else{
	  v_H[iz+Nz*(ix+Nx*iy)]=in[iy][ix][0][0]/(NNx*NNy*NNz)+_Complex_I*in[iy][ix][0][1]/(NNx*NNy*NNz);
	  }//*/
      }
    }
  } 
  
}


/**
 * 2. Second order gradients for grid functions
 * Using periodic condition 
 */

//template<class typename>
void gradient_real_2(double *f, double *fx, double *fy, double *fz, 
		     int Nx, int Ny, int Nz, double hx, double hy, double hz)
{
  int ix, iy, iz;
  // 1. peoridic second order gradient
    for(ix = 0; ix<Nx; ix++){
       for(iy = 0; iy<Ny; iy++){
          for(iz = 0; iz<Nz; iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+1+Nx*iy)] - f[iz+Nz*(Nx-2+Nx*iy)])/(2.0*hx);
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(1+Nx*iy)] - f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+1+Nx*iy)] - f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*(iy+1))] - f[iz+Nz*(ix+Nx*(Ny-2))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*1)] - f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*(iy+1))] - f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+1+Nz*(ix+Nx*iy)] - f[Nz-2+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[1+Nz*(ix+Nx*iy)] - f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+1+Nz*(ix+Nx*iy)] - f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/

}

void gradient_cplx_2(double complex *f, 
		     double complex *fx, double complex *fy, double complex *fz, 
		     int Nx, int Ny, int Nz, double hx, double hy, double hz)
{
  int ix, iy, iz;
  // 1. peoridic second order gradient
    for(ix = 0; ix<Nx; ix++){
       for(iy = 0; iy<Ny; iy++){
          for(iz = 0; iz<Nz; iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+1+Nx*iy)] - f[iz+Nz*(Nx-2+Nx*iy)])/(2.0*hx);
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(1+Nx*iy)] - f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+1+Nx*iy)] - f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*(iy+1))] - f[iz+Nz*(ix+Nx*(Ny-2))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*1)] - f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+Nz*(ix+Nx*(iy+1))] - f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+1+Nz*(ix+Nx*iy)] - f[Nz-2+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[1+Nz*(ix+Nx*iy)] - f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)] = 
		 (f[iz+1+Nz*(ix+Nx*iy)] - f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/

}


/**
 * 3. Fourth order
 * Using periodic condition 
 */

/*
template<class typename>
void grad_of_gridfun_4(typename *f, typename *fx, typename *fy, typename *fz, 
		      int Nx, int Ny, int Nz, double hx, double hy, double hz)
{
  int ix=0, iy=0, iz=0; 
  typename v0, v1, v2, v3;

  for(iy=0;iy<Ny;iy++){
    for(ix=0;ix<Nx;ix++){
      for(iz=0;iz<Nz;iz++){

	// x-derivatives
	if(ix==0){
	  f0 = wfs[iz+Nz*(Nx-3+Nx*iy)];
	  f1 = wfs[iz+Nz*(Nx-2+Nx*iy)];
	  f2 = wfs[iz+Nz*(1+Nx*(iy+Ny*occ))];
	  f3 = wfs[iz+Nz*(2+Nx*iy)]; 
	}
	else if(ix==1){
	  f0 = wfs[iz+Nz*(Nx-2+Nx*iy)];
	  f1 = wfs[iz+Nz*(0+Nx*iy)];
	  f2 = wfs[iz+Nz*(2+Nx*iy)];
	  f3 = wfs[iz+Nz*(3+Nx*iy)]; 
	}
	else if(ix==Nx-1){
	  f0 = wfs[iz+Nz*(Nx-3+Nx*iy)];
	  f1 = wfs[iz+Nz*(Nx-2+Nx*iy)];
	  f2 = wfs[iz+Nz*(1+Nx*iy)];
	  f3 = wfs[iz+Nz*(2+Nx*iy)]; 
	}
	else if(ix==Nx-2){
	  f0=wfs[iz+Nz*(Nx-4+Nx*iy)];
	  f1=wfs[iz+Nz*(Nx-3+Nx*iy)];
	  f2=wfs[iz+Nz*(Nx-1+Nx*iy)];
	  f3=wfs[iz+Nz*(1+Nx*iy)]; 
	}
	else{
	  f0=wfs[iz+Nz*(ix-2+Nx*iy)];
	  f1=wfs[iz+Nz*(ix-1+Nx*iy)];
	  f2=wfs[iz+Nz*(ix+1+Nx*iy)];
	  f3=wfs[iz+Nz*(ix+2+Nx*iy)];  
	}
	fx[iz+Nz*(ix+Nx*iy)] = (-f3+8*f2-8*f1+f0)/(12.0*hx); 

	// y-derivatives
	if(iy==0){
	  f0=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
	  f1=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
	  f2=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))];
	  f3=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))]; 
	}
	else if(iy==1){
	  f0=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
	  f1=wfs[iz+Nz*(ix+Nx*(0+Ny*occ))];
	  f2=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))];
	  f3=wfs[iz+Nz*(ix+Nx*(3+Ny*occ))]; 
	}
	else if(iy==Ny-1){
	  f0=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
	  f1=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
	  f2=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))];
	  f3=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))]; 
	}
	else if(iy==Ny-2){
	  f0=wfs[iz+Nz*(ix+Nx*(Ny-4+Ny*occ))];
	  f1=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
	  f2=wfs[iz+Nz*(ix+Nx*(Ny-1+Ny*occ))];
	  f3=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))]; 
	}
	else{
	  f0=wfs[iz+Nz*(ix+Nx*(iy-2+Ny*occ))];
	  f1=wfs[iz+Nz*(ix+Nx*(iy-1+Ny*occ))];
	  f2=wfs[iz+Nz*(ix+Nx*(iy+1+Ny*occ))];
	  f3=wfs[iz+Nz*(ix+Nx*(iy+2+Ny*occ))]; 
	}
	fy[iz+Nz*(ix+Nx*iy)] = (-f3+8*f2-8*f1+f0)/(12.0*hy);

	// z-dervatives
	if(iz==0){
	  f0=wfs[Nz-3+Nz*(ix+Nx*iy)];
	  f1=wfs[Nz-2+Nz*(ix+Nx*iy)];
	  f2=wfs[1+Nz*(ix+Nx*iy)];
	  f3=wfs[2+Nz*(ix+Nx*iy)]; 
	}
	else if(iz==1){
	  f0=wfs[Nz-2+Nz*(ix+Nx*iy)];
	  f1=wfs[0+Nz*(ix+Nx*iy)];
	  f2=wfs[2+Nz*(ix+Nx*iy)];
	  f3=wfs[3+Nz*(ix+Nx*iy)]; 
	}
	else if(iz==Nz-1){
	  f0=wfs[Nz-3+Nz*(ix+Nx*iy)];
	  f1=wfs[Nz-2+Nz*(ix+Nx*iy)];
	  f2=wfs[1+Nz*(ix+Nx*iy)];
	  f3=wfs[2+Nz*(ix+Nx*iy)]; 
	}
	else if(iz==Nz-2){
	  f0=wfs[Nz-4+Nz*(ix+Nx*iy)];
	  f1=wfs[Nz-3+Nz*(ix+Nx*iy)];
	  f2=wfs[Nz-1+Nz*(ix+Nx*iy)];
	  f3=wfs[1+Nz*(ix+Nx*iy)]; 
	}
	else{
	  f0=wfs[iz-2+Nz*(ix+Nx*iy)];
	  f1=wfs[iz-1+Nz*(ix+Nx*iy)];
	  f2=wfs[iz+1+Nz*(ix+Nx*iy)];
	  f3=wfs[iz+2+Nz*(ix+Nx*iy)]; 
	}
	fz[iz+Nz*(ix+Nx*iy)] = (-f3+8*f2-8*f1+f0)/(12.0*hz);

      }    // ix 
    }   // iz
  }  // iy

} // end of gradient_grid_function

*/


void gradient_of_wfs(double *wfs, double *gradient_wfs, int st_num, 
		     int Nx, int Ny, int Nz, 
		     double hx, double hy, double hz)
{
  int ix=0, iy=0, iz=0;
  int occ=0;

  double f0, f1, f2, f3;

  //* Fourth order
  for(ix=0;ix<Nx;ix++){
     for(iy=0;iy<Ny;iy++){
        for(iz=0;iz<Nz;iz++){
           //Using periodic condition 
           for(occ=0;occ<st_num;occ++){
              if(ix==0){
                f0=wfs[iz+Nz*(Nx-3+Nx*(iy+Ny*occ))];
                f1=wfs[iz+Nz*(Nx-2+Nx*(iy+Ny*occ))];
                f2=wfs[iz+Nz*(1+Nx*(iy+Ny*occ))];
                f3=wfs[iz+Nz*(2+Nx*(iy+Ny*occ))]; 
              }
              else if(ix==1){
                f0=wfs[iz+Nz*(Nx-2+Nx*(iy+Ny*occ))];
                f1=wfs[iz+Nz*(0+Nx*(iy+Ny*occ))];
                f2=wfs[iz+Nz*(2+Nx*(iy+Ny*occ))];
                f3=wfs[iz+Nz*(3+Nx*(iy+Ny*occ))]; 
              }
              else if(ix==Nx-1){
                f0=wfs[iz+Nz*(Nx-3+Nx*(iy+Ny*occ))];
                f1=wfs[iz+Nz*(Nx-2+Nx*(iy+Ny*occ))];
                f2=wfs[iz+Nz*(1+Nx*(iy+Ny*occ))];
                f3=wfs[iz+Nz*(2+Nx*(iy+Ny*occ))]; 
              }
              else if(ix==Nx-2){
                f0=wfs[iz+Nz*(Nx-4+Nx*(iy+Ny*occ))];
                f1=wfs[iz+Nz*(Nx-3+Nx*(iy+Ny*occ))];
                f2=wfs[iz+Nz*(Nx-1+Nx*(iy+Ny*occ))];
                f3=wfs[iz+Nz*(1+Nx*(iy+Ny*occ))]; 
              }
              else{
                f0=wfs[iz+Nz*(ix-2+Nx*(iy+Ny*occ))];
                f1=wfs[iz+Nz*(ix-1+Nx*(iy+Ny*occ))];
                f2=wfs[iz+Nz*(ix+1+Nx*(iy+Ny*occ))];
                f3=wfs[iz+Nz*(ix+2+Nx*(iy+Ny*occ))];  
              }
              gradient_wfs[0+3*(iz+Nz*(ix+Nx*(iy+Ny*occ)))]=(-f3+8*f2-8*f1+f0)/(12.0*hx); 


              if(iy==0){
                f0=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
                f1=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
                f2=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))];
                f3=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))]; 
              }
              else if(iy==1){
                f0=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
                f1=wfs[iz+Nz*(ix+Nx*(0+Ny*occ))];
                f2=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))];
                f3=wfs[iz+Nz*(ix+Nx*(3+Ny*occ))]; 
              }
              else if(iy==Ny-1){
                f0=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
                f1=wfs[iz+Nz*(ix+Nx*(Ny-2+Ny*occ))];
                f2=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))];
                f3=wfs[iz+Nz*(ix+Nx*(2+Ny*occ))]; 
              }
              else if(iy==Ny-2){
                f0=wfs[iz+Nz*(ix+Nx*(Ny-4+Ny*occ))];
                f1=wfs[iz+Nz*(ix+Nx*(Ny-3+Ny*occ))];
                f2=wfs[iz+Nz*(ix+Nx*(Ny-1+Ny*occ))];
                f3=wfs[iz+Nz*(ix+Nx*(1+Ny*occ))]; 
              }
              else{
                f0=wfs[iz+Nz*(ix+Nx*(iy-2+Ny*occ))];
                f1=wfs[iz+Nz*(ix+Nx*(iy-1+Ny*occ))];
                f2=wfs[iz+Nz*(ix+Nx*(iy+1+Ny*occ))];
                f3=wfs[iz+Nz*(ix+Nx*(iy+2+Ny*occ))]; 
              }
              gradient_wfs[1+3*(iz+Nz*(ix+Nx*(iy+Ny*occ)))]=(-f3+8*f2-8*f1+f0)/(12.0*hy);

              if(iz==0){
                f0=wfs[Nz-3+Nz*(ix+Nx*(iy+Ny*occ))];
                f1=wfs[Nz-2+Nz*(ix+Nx*(iy+Ny*occ))];
                f2=wfs[1+Nz*(ix+Nx*(iy+Ny*occ))];
                f3=wfs[2+Nz*(ix+Nx*(iy+Ny*occ))]; 
              }
              else if(iz==1){
                f0=wfs[Nz-2+Nz*(ix+Nx*(iy+Ny*occ))];
                f1=wfs[0+Nz*(ix+Nx*(iy+Ny*occ))];
                f2=wfs[2+Nz*(ix+Nx*(iy+Ny*occ))];
                f3=wfs[3+Nz*(ix+Nx*(iy+Ny*occ))]; 
              }
              else if(iz==Nz-1){
                f0=wfs[Nz-3+Nz*(ix+Nx*(iy+Ny*occ))];
                f1=wfs[Nz-2+Nz*(ix+Nx*(iy+Ny*occ))];
                f2=wfs[1+Nz*(ix+Nx*(iy+Ny*occ))];
                f3=wfs[2+Nz*(ix+Nx*(iy+Ny*occ))]; 
              }
              else if(iz==Nz-2){
                f0=wfs[Nz-4+Nz*(ix+Nx*(iy+Ny*occ))];
                f1=wfs[Nz-3+Nz*(ix+Nx*(iy+Ny*occ))];
                f2=wfs[Nz-1+Nz*(ix+Nx*(iy+Ny*occ))];
                f3=wfs[1+Nz*(ix+Nx*(iy+Ny*occ))]; 
              }
              else{
                f0=wfs[iz-2+Nz*(ix+Nx*(iy+Ny*occ))];
                f1=wfs[iz-1+Nz*(ix+Nx*(iy+Ny*occ))];
                f2=wfs[iz+1+Nz*(ix+Nx*(iy+Ny*occ))];
                f3=wfs[iz+2+Nz*(ix+Nx*(iy+Ny*occ))]; 
              }
              gradient_wfs[2+3*(iz+Nz*(ix+Nx*(iy+Ny*occ)))]=(-f3+8*f2-8*f1+f0)/(12.0*hz);
           }
        }
     }
  }//*/
  
  //*/

  return;
}


void gradient_hessian_eta_xc(double complex *eta_xc,
			     double complex *gradient_eta_xc, 
			     double complex *hessian_eta_xc,  
			     int Nx, int Ny, int Nz, 
			     double hx, double hy, double hz)
{

  int ix=0, iy=0, iz=0; 
  double complex v0, v1, v2, v3;
  int icase;
  
  icase=2;
  /********************
  icase=1: non periodic
  icase=2: periodic
  *********************/
  if(icase==1){
  }
  else if(icase==2){
    /*****************************
    calculate gradient of eta_xc
    ******************************/
    //* Second order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){
               gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(ix+1+Nx*iy)]-eta_xc[iz+Nz*(Nx-2+Nx*iy)])/(2.0*hx);
             }
             else if(ix==Nx-1){
               gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(1+Nx*iy)]-eta_xc[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
             else{
               gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(ix+1+Nx*iy)]-eta_xc[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(ix+Nx*(iy+1))]-eta_xc[iz+Nz*(ix+Nx*(Ny-2))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(ix+Nx*1)]-eta_xc[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }
             else{
               gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+Nz*(ix+Nx*(iy+1))]-eta_xc[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+1+Nz*(ix+Nx*iy)]-eta_xc[Nz-2+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else if(iz==Nz-1){
               gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[1+Nz*(ix+Nx*iy)]-eta_xc[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else{
               gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*iy))]=(eta_xc[iz+1+Nz*(ix+Nx*iy)]-eta_xc[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/
    
    /********************************
    Calculate hessian_eta_xc
    ********************************/  
    //* Second order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){ 
               hessian_eta_xc[0+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[0+3*(iz+Nz*(Nx-2+Nx*iy))])/(2.0*hx); 
               hessian_eta_xc[0+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[1+3*(iz+Nz*(Nx-2+Nx*iy))])/(2.0*hx); 
               hessian_eta_xc[0+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[2+3*(iz+Nz*(Nx-2+Nx*iy))])/(2.0*hx); 
             }
             else if(ix==Nx-1){
               hessian_eta_xc[0+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(1+Nx*iy))]-gradient_eta_xc[0+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
               hessian_eta_xc[0+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(1+Nx*iy))]-gradient_eta_xc[1+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
               hessian_eta_xc[0+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(1+Nx*iy))]-gradient_eta_xc[2+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
             }
             else{
               hessian_eta_xc[0+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[0+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
               hessian_eta_xc[0+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[1+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
               hessian_eta_xc[0+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(ix+1+Nx*iy))]-gradient_eta_xc[2+3*(iz+Nz*(ix-1+Nx*iy))])/(2.0*hx);
             }
           
             if(iy==0){
               hessian_eta_xc[1+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*(Ny-2)))])/(2.0*hy);
               hessian_eta_xc[1+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*(Ny-2)))])/(2.0*hy);
               hessian_eta_xc[1+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*(Ny-2)))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               hessian_eta_xc[1+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*1))]-gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
               hessian_eta_xc[1+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*1))]-gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
               hessian_eta_xc[1+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*1))]-gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
             }
             else{
               hessian_eta_xc[1+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[0+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
               hessian_eta_xc[1+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[1+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
               hessian_eta_xc[1+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*(iy+1)))]-gradient_eta_xc[2+3*(iz+Nz*(ix+Nx*(iy-1)))])/(2.0*hy);
             }

             if(iz==0){
               hessian_eta_xc[2+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[0+3*(Nz-2+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[1+3*(Nz-2+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[2+3*(Nz-2+Nz*(ix+Nx*iy))])/(2.0*hz);
             }
             else if(iz==Nz-1){
               hessian_eta_xc[2+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(1+Nz*(ix+Nx*iy))]-gradient_eta_xc[0+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(1+Nz*(ix+Nx*iy))]-gradient_eta_xc[1+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(1+Nz*(ix+Nx*iy))]-gradient_eta_xc[2+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
             }
             else{
               hessian_eta_xc[2+3*(0+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[0+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[0+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(1+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[1+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[1+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
               hessian_eta_xc[2+3*(2+3*(iz+Nz*(ix+Nx*iy)))]=(gradient_eta_xc[2+3*(iz+1+Nz*(ix+Nx*iy))]-gradient_eta_xc[2+3*(iz-1+Nz*(ix+Nx*iy))])/(2.0*hz);
             }
          }
       }
    }//*/
  }
  else{
    printf(" Wrong choice of gradient_hessian of eta_xc \n"); getchar();
  }
  return;
}


void gradient_of_complex_function(double complex *f, 
				  double complex *fx, 
				  double complex *fy, 
				  double complex *fz, 
				  int Nx, int Ny, int Nz, 
				  double hx, double hy, double hz)
{
  int ix=0, iy=0, iz=0; 
  double complex v0, v1, v2, v3;

  int icase=0;
  /**********************
  icase= 1, non periodic
  icase= 2, periodic
  ***********************/
  icase=2;
  
  if(icase==1){
    //* First order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix+Nx*iy)])/hx;
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/hx;
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*iy)])/hy;
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix+Nx*(iy-1))])/hy;
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix+Nx*iy)])/hz;
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/hz;
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/
  }
  else if(icase==2){
    //* Second order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(Nx-2+Nx*iy)])/(2.0*hx);
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(Ny-2))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*1)]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[Nz-2+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)]=(f[1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/
  }
  else{
    printf(" Wrong choice of gradient of complex function \n "); getchar();
  }
  return;
}

// for derivatives of rho_0
void gradient_of_real_function(double  *f, double  *fx, double  *fy, double  *fz, 
			       int Nx, int Ny, int Nz, 
			       double hx, double hy, double hz)
{
  int ix=0, iy=0, iz=0; 
  double v0, v1, v2, v3;

  int icase;
  /**********************
  icase= 1, non periodic
  icase= 2, periodic
  ***********************/
  icase=2;
  
  if(icase==1){
    //* First order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix+Nx*iy)])/hx;
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/hx;
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*iy)])/hy;
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix+Nx*(iy-1))])/hy;
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz+Nz*(ix+Nx*iy)])/hz;
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/hz;
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/
  }
  else if(icase==2){
    //* Second order
    for(ix=0;ix<Nx;ix++){
       for(iy=0;iy<Ny;iy++){
          for(iz=0;iz<Nz;iz++){
             if(ix==0){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(Nx-2+Nx*iy)])/(2.0*hx);
             }
             else if(ix==Nx-1){
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
             else{
               fx[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+1+Nx*iy)]-f[iz+Nz*(ix-1+Nx*iy)])/(2.0*hx);
             }
           
             if(iy==0){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(Ny-2))])/(2.0*hy);
             }
             else if(iy==Ny-1){
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*1)]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }
             else{
               fy[iz+Nz*(ix+Nx*iy)]=(f[iz+Nz*(ix+Nx*(iy+1))]-f[iz+Nz*(ix+Nx*(iy-1))])/(2.0*hy);
             }

             if(iz==0){
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[Nz-2+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else if(iz==Nz-1){
               fz[iz+Nz*(ix+Nx*iy)]=(f[1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
             else{
               fz[iz+Nz*(ix+Nx*iy)]=(f[iz+1+Nz*(ix+Nx*iy)]-f[iz-1+Nz*(ix+Nx*iy)])/(2.0*hz);
             }
          }
       }
    }//*/
  }
  else{
    printf(" Wrong choice of gradient of real function \n "); getchar();
  }
  return;
}


// interpolating subroutine 1
void mesh_2_grid(double complex *AEM_1, double complex *AEM_2,  
		 double complex *AEM_3,   
		 double *x1, double *y1, double *z1, int Nx, int Ny, int Nz,
		 double *sampleX, double *sampleY, double *sampleZ, int N_sample, 
		 double complex *A_EM, 
		 double *Node, int *Element, int *Element_Edge, 
		 double *Length_of_Edge, 
		 int *PointInTet)
{
  int isample, ix,iy,iz;
  int e, ie, i1,i2,ii,jj,mm,nn;

  double x0,   y0,  z0;
  double xx0, yy0, zz0;
  double xe[4],ye[4],ze[4];
  double ae[4],be[4],ce[4],de[4];
  double Ve, Lijmn[4];
  double aem1, aem2, aem3;

  int Edge_end[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};  // in local index

  for(isample=0;isample<N_sample;isample++){
    x0=sampleX[isample]; y0=sampleY[isample]; z0=sampleZ[isample];

    for(ix=0;ix<Nx;ix++){
      for(iy=0;iy<Ny;iy++){
	for(iz=0;iz<Nz;iz++){ 
	  xx0 = x0 + x1[ix]; yy0 = y0 + y1[iy]; zz0 = z0 + z1[iz];
	  e = PointInTet[iz+Nz*(ix+Nx*(iy+Ny*isample))]-1; 
	  ii=Element[0+4*e]-1; 
	  jj=Element[1+4*e]-1; 
	  mm=Element[2+4*e]-1; 
	  nn=Element[3+4*e]-1;  

	  xe[0]=Node[0+3*ii]; xe[1]=Node[0+3*jj]; xe[2]=Node[0+3*mm]; xe[3]=Node[0+3*nn];
	  ye[0]=Node[1+3*ii]; ye[1]=Node[1+3*jj]; ye[2]=Node[1+3*mm]; ye[3]=Node[1+3*nn];
	  ze[0]=Node[2+3*ii]; ze[1]=Node[2+3*jj]; ze[2]=Node[2+3*mm]; ze[3]=Node[2+3*nn];

	  ae[0]= xe[1]*(ye[2]*ze[3]-ze[2]*ye[3])-xe[2]*(ye[1]*ze[3]-ze[1]*ye[3])+xe[3]*(ye[1]*ze[2]-ze[1]*ye[2]);
	  ae[1]=-xe[0]*(ye[2]*ze[3]-ze[2]*ye[3])+xe[2]*(ye[0]*ze[3]-ze[0]*ye[3])-xe[3]*(ye[0]*ze[2]-ze[0]*ye[2]);
	  ae[2]= xe[0]*(ye[1]*ze[3]-ze[1]*ye[3])-xe[1]*(ye[0]*ze[3]-ze[0]*ye[3])+xe[3]*(ye[0]*ze[1]-ze[0]*ye[1]);
	  ae[3]=-xe[0]*(ye[1]*ze[2]-ze[1]*ye[2])+xe[1]*(ye[0]*ze[2]-ze[0]*ye[2])-xe[2]*(ye[0]*ze[1]-ze[0]*ye[1]);

	  be[0]= (ye[2]*ze[3]-ze[2]*ye[3])-(ye[1]*ze[3]-ze[1]*ye[3])+(ye[1]*ze[2]-ze[1]*ye[2]);
	  be[1]=-(ye[2]*ze[3]-ze[2]*ye[3])+(ye[0]*ze[3]-ze[0]*ye[3])-(ye[0]*ze[2]-ze[0]*ye[2]);
	  be[2]= (ye[1]*ze[3]-ze[1]*ye[3])-(ye[0]*ze[3]-ze[0]*ye[3])+(ye[0]*ze[1]-ze[0]*ye[1]);
	  be[3]=-(ye[1]*ze[2]-ze[1]*ye[2])+(ye[0]*ze[2]-ze[0]*ye[2])-(ye[0]*ze[1]-ze[0]*ye[1]);
	  be[0]=-be[0]; be[1]=-be[1]; be[2]=-be[2]; be[3]=-be[3];

	  ce[0]= (xe[2]*ze[3]-ze[2]*xe[3])-(xe[1]*ze[3]-ze[1]*xe[3])+(xe[1]*ze[2]-ze[1]*xe[2]);
	  ce[1]=-(xe[2]*ze[3]-ze[2]*xe[3])+(xe[0]*ze[3]-ze[0]*xe[3])-(xe[0]*ze[2]-ze[0]*xe[2]);
	  ce[2]= (xe[1]*ze[3]-ze[1]*xe[3])-(xe[0]*ze[3]-ze[0]*xe[3])+(xe[0]*ze[1]-ze[0]*xe[1]);
	  ce[3]=-(xe[1]*ze[2]-ze[1]*xe[2])+(xe[0]*ze[2]-ze[0]*xe[2])-(xe[0]*ze[1]-ze[0]*xe[1]);

	  de[0]= (xe[2]*ye[3]-ye[2]*xe[3])-(xe[1]*ye[3]-ye[1]*xe[3])+(xe[1]*ye[2]-ye[1]*xe[2]);
	  de[1]=-(xe[2]*ye[3]-ye[2]*xe[3])+(xe[0]*ye[3]-ye[0]*xe[3])-(xe[0]*ye[2]-ye[0]*xe[2]);
	  de[2]= (xe[1]*ye[3]-ye[1]*xe[3])-(xe[0]*ye[3]-ye[0]*xe[3])+(xe[0]*ye[1]-ye[0]*xe[1]);
	  de[3]=-(xe[1]*ye[2]-ye[1]*xe[2])+(xe[0]*ye[2]-ye[0]*xe[2])-(xe[0]*ye[1]-ye[0]*xe[1]);
	  de[0]=-de[0]; de[1]=-de[1]; de[2]=-de[2]; de[3]=-de[3];

	  Ve=(ae[0]+ae[1]+ae[2]+ae[3])/6.0;    
	  for(ii=0;ii<4;ii++){ 
	    ae[ii]=ae[ii]/(6.0*Ve); 
	    be[ii]=be[ii]/(6.0*Ve);
	    ce[ii]=ce[ii]/(6.0*Ve); 
	    de[ii]=de[ii]/(6.0*Ve);
	  }
 
	  Lijmn[0]=ae[0]+be[0]*xx0+ce[0]*yy0+de[0]*zz0;
	  Lijmn[1]=ae[1]+be[1]*xx0+ce[1]*yy0+de[1]*zz0;
	  Lijmn[2]=ae[2]+be[2]*xx0+ce[2]*yy0+de[2]*zz0;
	  Lijmn[3]=ae[3]+be[3]*xx0+ce[3]*yy0+de[3]*zz0; 

	  aem1=0; aem2=0; aem3=0;
	  for(ie=0;ie<6;ie++){ 
	    i1=Edge_end[ie][0]; i2=Edge_end[ie][1];
	    aem1=aem1+A_EM[Element_Edge[ie+6*e]-1]*(Lijmn[i1]*be[i2]-Lijmn[i2]*be[i1])*Length_of_Edge[Element_Edge[ie+6*e]-1];
	    aem2=aem2+A_EM[Element_Edge[ie+6*e]-1]*(Lijmn[i1]*ce[i2]-Lijmn[i2]*ce[i1])*Length_of_Edge[Element_Edge[ie+6*e]-1];
	    aem3=aem3+A_EM[Element_Edge[ie+6*e]-1]*(Lijmn[i1]*de[i2]-Lijmn[i2]*de[i1])*Length_of_Edge[Element_Edge[ie+6*e]-1];
	  } 
	  AEM_1[iz+Nz*(ix+Nx*(iy+Ny*isample))]=aem1;
	  AEM_2[iz+Nz*(ix+Nx*(iy+Ny*isample))]=aem2;
	  AEM_3[iz+Nz*(ix+Nx*(iy+Ny*isample))]=aem3;  

	}    //iz
      }    //iy
    }    //ix

  }    // isample

}

// interpolating subroutine 2
void grid_2_mesh(double complex *del_j ,    // out put FEM function
		 double *Node, int *Element, int *Element_Edge, 
		 int N_edge, double *Length_of_Edge,
		 double complex *del_j_1,   // input variables
		 double complex *del_j_2, 
		 double complex *del_j_3,
		 double *x1, double *y1, double *z1, int Nx, int Ny, int Nz,
		 double *sampleX, double *sampleY, double *sampleZ, int N_sample, 
		 int *NodeInSamples)
{
  
}

